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(Received ) 

We describe a quantum algorithm to compute the density of states and thermal equilib- 
rium properties of quantum many-body systems. We present results obtained by running 
this algorithm on a software implementation of a 21-qubit quantum computer for the case 
of an antiferromagnetic Heisenberg model on triangular lattices of different size. 

§1. Introduction 



Q> ' Recent theoretical work has shown that a Quantum Computer (QC) has the 

. potential of solving certain computationally hard problems such as factoring integers 

Oh! and searching databases much faster than a conventional computer. " The idea 

that a QC might be more powerful than an ordinary computer is based on the notion 
that a quantum system can be in any superposition of states and that interference of 
^ ' these states allows exponentially many computations to be done in parallel. This 

hypothetical power of a QC might be used to solve other difficult problems as well, 
such as for example the calculation of the physical properties of quantum many- 
^ \ body systems. " As a matter of fact, part of Feynman's original motivation to 

^ ■ consider QC's was that they might be used as a vehicle to perform exact simulations 

of quantum mechanical phenomena. For future applications it is clearly of interest 
to address the question how to program a QC such that it performs a simulation of 
specific physical systems. 

In this paper we describe a quantum algorithm (QA) to compute the equilib- 
rium properties of quantum systems by making a few statistically uncorrelated runs, 
starting from random initial states. Exploiting the intrinsic parallellism of the hy- 
pothetical QC this QA executes in polynomial time. En route this QA computes 
the density of states (DOS) from which all the eigenvalues of the model Hamiltonian 
may be determined. We test our QA on a software implemention of a 21-qubit QC 
by explicit calculations for an antiferromagnetic Heisenberg model on a triangular 
lattice. 



E-mail: deraedt@phys.rug.nl 



2 



H. De Raedt et al. 



§2. Theory 

Consider the problem of computing the physical properties of a quantum system, 
described by a Hamiltonian H, in thermal equilibrium at a temperature T. In the 
canonical ensemble this equilibrium state is characterised by the partition function 
Z = Trexp{—PH) where /3 = l/ksT and ks is Boltzmann's constant (we put ks = l 
and ^ = 1 in the rest of this paper). The dimension of the Hilbert space of physical 
states will be denoted hy D. If all the eigenvalues {Ei ;{ = !,..., D} of H are known, 
we can make use of the fact that Z = J2iLi 6xp(— /^-Bi) to reduce the computation 
of the physical properties to a classical statistical mechanical problem which can be 
solved by standard probabilistic methods. For a non-trivial quantum many-body 
system the determination of the eigenvalues is a difficult computational problem 
itself, in practice as difficult as the calculation of Z. In view of this we will from now 
on assume that the eigenvalues of H are not known. 

The DOS 

1 r+oo 

V{e) = E '^(^ - ^0 = 2^ '^e"'*'' dt, (1) 

determines the equilibrium state of the system. Indeed, once P(e) is known Z is 
easy to calculate: 



Z 



/ + 00 
e-^'V{€) de. (2) 
-oo 



Integral (2) exists whenever the spectrum of H has a lowerbound, i.e. P(e) = for 
— oo < e < minjE^j. Note that 'D{e) is a real-valued function. 

With suitable (time-dependent) modifications of H, the partition function plays 
the role of a generating function from which all physical quantities of interest can be 
obtained. Physical quantities such as the energy and specific heat are given by 

1 r+oo 

E = {H) = - ee-^^V{e)de, (3) 

^ J —oo 



and 



C = /3\{H^) - {Hf) = e^e-^^V{e) de - E^'^ , (4) 



respectively. 

From (1) it is clear that the computation of T'(e) consists of two parts: 1) 
Compute the trace of e~**^ for many values of t and 2) perform a (Fast) Fourier 
Transform of this data. It is known how to carry out part 2 on a QC " so we 
focus on part 1. 

The QA described below computes e~^^^\ip) for any ^ in O (logD) operations. 
We now argue that in practice it will usually be sufficient to determine a small 
fraction (?» O (logD)) of the D diagonal matrix elements of e~**^. Instead of com- 
puting diagonal matrix elements with respect to a chosen complete set of basis states 
{(pn, n = 1, . . . , D}, we generate random numbers {a„; n = 1, . . . , D} and construct 
the new state 

D 

|^) = E«n|'/'n). (5) 
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The corresponding diagonal element of the time-evolution operator reads 

(#|e-^*^#) = J2 a*^aMe-''''<Pm). (6) 

n,m=l 

If we now generate the a^'s such that a*a^ = 6n,m denotes the average of x over 
statistically independent realizations) then 

D 

{^\e-iiH^) = ^(,/)„|e-^*^</.„) = TVe-'*^. (7) 

n=l 

In other words, the trace of the time-evolution operator can be estimated by random 
sampling of the states |^). 

A QC computes \e~^*^^) just as easily as |e~**^i;/)m) and in practice there 
would be no need to have a random state generator: Switching the QC off and 
on will put the QC in some random initial state. However to compute (^|e~**^^) 
we would have to store the initial state (|<P)) in the QC and project |e~'*^#) 
onto it, a rather complicated procedure. Instead it is more effective to apply to 
{(pi) a random sequence of Controlled-NOT operations to construct e.g. an en- 
tangled random state \0) = D~-^/^(zb^i it 02 . • • ± (po)-^^^'^'^^ We then calculate 
|e-'*^/2^) = En=i&n(i/2) and use 

D 

(^|e-^*^^) = (e^*^/2^|e-^*^/2#) = ^ \bn{t/2)f (8) 

n=l 

to obtain the diagonal matrix element. Each of these steps executes very efficiently 
on a QC. 

Remains the question how many samples S are needed to compute the energy 
and specific heat to high accuracy. According to the central limit theorem the statis- 
tical error on the results vanishes as 1/^/S. However, as we demonstrate below, the 
application of our QA to a highly non-trivial quantum many-body system provides 
strong evidence that this error also decreases with the system size. For systems 
of 15 qubits or more, we find that taking S" = 20 samples already gives very accu- 
rate results. Our experimental finding that the statistical error on (^|e-»*^^) for 
randomly chosen |#) decreases with D gives an extra boost to the efficiency of the 
QA. 

§3. Soft Quantum Computer and Quantum Algorithm 

The method described above has been tested on our Soft Quantum Computer 
(SQC). The SQC used to compute the results presented in the present paper is a 
hard-coded version, derived from of a more versatile SQC discussed elsewhere. 
Our SQC solves the time-dependent Schrodinger equation (TDSE) 

i^\m) = H\m), (9) 
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for a quantum many-body system described by the spin-1/2 Hamiltonian 

i,j=l a=x,y,z i=l a=x,y,z 

where the first sum runs over all pairs P of spins (qubits), Sf^ (a = x, y, z) denotes the 
a-th component of the spin-1/2 operator representing the i-th spin, J^j determines 
the strength of the interaction between the spins at sites i and j, and hf is the 
(local magnetic) field acting on the i-th spin. The number of qubits is L and D = 
2^. Hamiltonian (10) is sufficiently general to capture the salient features of most 
physical models of QC's (our SQC also deals with time-dependent external fields). 

According to (9) the QC will evolve in time through the D x D unitary transfor- 
mation U{t) = e"**^. We now describe the QA that compTitcs U{t)\'P) for arbitrary 
1^). Using the semi-group property of U{t) to write U{t) = U{t)"^ where t = mr, 
the main step is to replace U{t) by a symmetrized product-formula approxima- 
tion. For the case at hand it is expedient to take 

U{t) « U{t) = e-irH,/^e-'^"y/^e-'^"-e-'''^y/^e-'^"^/^, (11) 

where 

L L 

Ha = -J2 -^^J^i^J - E htSr ; a = x,y,z. (12) 

i,j=l 1=1 

Evidently U{t) is unitary and hence the algorithm to solve the TDSE is uncondi- 
tionally stable. 

As basis states {|<^n)} we take the direct product of the eigenvectors of the Sf 
(i.e. spin-up and spin-down |ii)). In this basis, e~^'^^'/^ changes the input 
state by altering the phase of each of the basis vectors. As is a sum of pair 
interactions it is trivial to rewrite this operation as a direct product of 4x4 diagonal 
matrices (containing the interaction-controlled phase shifts) and 4x4 unit matrices. 
Still working in the same representation, the action of e~*'^^^/^ can be written in a 
similar manner but the matrices that contain the interaction-controlled phase-shift 
have to be replaced by non-diagonal matrices. Although this docs not present a real 
problem it is more efficient and systematic to proceed as follows. Let us denote by 
X{Y) the rotation by 7r/2 of each spin about the a;(y)-axis. As 

^-iTHy/2 ^ xX^e-'^^y/^XX^ = Xe-'^^'^'^X\ (13) 

it is clear that the action of e~*^^^/^ can be computed by applying to each qubit, 
the inverse of X followed by an interaction-controlled phase-shift and X itself. The 
prime in (13) indicates that Jf^ and /t| in Hz have to be replaced by J^^ and h\ 
respectively. A similar procedure is used to compute the action of e"*"^^^ . 

Our SQC carries out O (^P 2^^ operations to perform the transformation e^*"^^^/^ 
but a QC operates on all qubits simultaneously and would therefore only need O (P) 
operations. The operation counts for e"*'^-'^^ (or e~^'^^y ) are O ^(P -|- 2)2^^ and 

O {P + 2) for the SQC and QC respectively. On a QC the total operation count per 
time-step is 0(3P-|-4). 
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Fig. 1. Specific heat per site as a function of 
the temperature. 



Fig. 2. Standard deviation on the specific heat 
per site. 



§4. Application 



The QA described above has been tested on our SQC by simulating the an- 

tiferromagnctic spin 1/2 Heisenberg model with J = — 1 on triangular lattices of 
L = 6, 10, 15, 21 sites, subject to free boundary conditions. The ground-state prop- 
erties of this model can be computed by standard sparse-matrix techniques, see e.g. 
Ref. 20). The low temperature properties of this model are difficult to compute by 
conventional Quantum Monte Carlo (QMC) methods. The presence of frustrated 
interactions leads to the sign problem that is very often encountered in QMC 

work. 22). 23) 

In Fig. 1 we present some SQC results for the specific heat per site E/L as 

a function of the temperature. The number of samples S* = 20 in all cases. Also 
shown is data obtained by exact diagonalization of H iov L = 6, 10. *) The SQC and 
exact results differ significantly for temperatures where the specific heat exhibits 
a sharp peak. This is related to the presence of a gap in the low-energy part of 
the spectrum and the random choice of the |cp)'s. In this low-temperature regime 
where only a few of the lowest eigenvalues contribute, random fluctuations can have 
a large effect. However this is not really a problem: Knowing the DOS it is not 
difficult to determine the precise values of these few eigenvalues and compute C 
directly, without invoking (4). In Fig. 2 we show results for the standard deviation 
(SD) on C/L, calculated from the same data. The most remarkable feature of the 
SD is its dependence on the system size: The larger the system, the smaller the 
SD. Unfortunately we cannot yet offer a theoretical basis for this observed decrease. 
At very low temperature the SD on C/L goes to zero because only one eigenstate 
(i.e. the ground state) effectively contributes. The large values of the L = 6, 10 
low temperature SD data reflect the fact that in this regime the SQC and exact 
results differ considerably and indicate that for these system sizes more than S = 20 
samples are necessary to obtain accurate results. On the other hand, except for test 
purposes, we wouldn't use a QC to simulate a 10-site system because it can easily 



*^ The calculation of all 32768 (2097152) eigenvalues of the L = 15 (L = 21) system by standard 
linear algebra methods requires tremendous computational resources 
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be solved exactly on an ordinary computer. 



§5. Summary 



We have described a quantum algorithm to determine the distribution of eigen- 
values of a quantum many-body system in polynomial time. Prom these data thermal 
equilibrium properties of the system can be computed directly. The approach has 
been illustrated by numerical calculations on a software emulator of a physical model 
of a quantum computer. Excellent results have been obtained, suggesting a new route 
for simulating experiments on quantum systems on a quantum computer. However, 
implicit in the formulation of this physical model of the quantum computer is the 
assumption that each physical spin represents one qTibit. If this were not the case, 
the quantum computer will operate with much less efficiency. 
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